2.2.1 Métadonnées

Le dossier des limites administratives du Benin fourni par HDX (cliquer ici) est accessible au format shapefile et comporte différentes couches correspondant aux différents niveaux administratifs. Sa datation laisse augurer qu’il s’agit de la situation la plus récente: les fichiers datent de juin 2021 mais leur nom comporte la séquence “20190816” qui doit plutôt correspondre au 16 Août 2019. En pratique, les limites sont celles du recensement de 2013.

N.B. nous avons toutefois découvert une erreur importante dans le code de la commune d’Adjara qui est noté BJ0101 alors que le code exact est BJ1001. Nous avons également constaté une variabilité des noms de vraiables désignant les codes dans la base OCHA. C’est pourquoi nous avons préféré harmoniser les fichiers.

adm0<-st_read("data/admin/ben_adm0.shp",quiet = TRUE)
  adm1<-st_read("data/admin/ben_adm1.shp",quiet = TRUE)
  adm2<-st_read("data/admin/ben_adm2.dbf", quiet = TRUE)

2.2.2 Exemple d’application

Pour vérifier la qualité des fichiers nous effectuons un calage sur un fonds de carte OpenStreetMap à l’aide du package leaflet :

map <- leaflet() %>% 
              addTiles() %>%
              setView(lat = 8, lng=2.1, zoom = 6) %>%
                addPolygons(data = adm2,
                            label = ~ADM2_NAME,
                            highlightOptions = TRUE,
                             fill = TRUE,
                             color = "orange",
                             weight = 3 ) %>%
               addPolygons(data = adm1,
                          fill = FALSE,
                          color = "red",
                          weight = 2) %>%
               addPolygons(data = adm0,
                          fill = FALSE,
                          color = "black",
                          weight = 1) 
    
  
  map

2.3.1 Métadonnées

Un dossier du HDX accessible ici fournit les populations par âge et par sexe en 2019 au niveau des communes, départements ou pays entier. Il s’agit naturellement d’estimations mais les données n’en sont pas moins très utile pour toute analyse travaillant sur les situations présentes. Les classes d’âges sont suffisamment détaillées (tranches de 5 ans de 0 à 80 ans) pour procéder à des analyses démographiques intéressantes.

2.3.2 Exemple d’application

A titre de vérification, on réalise une analyse factorielle des correspondances AFC sur la pyramide des âges des 77 communes

don<-read.table("data/pop2019/ben_adm2_pop.csv", sep=";",header=T, encoding = "UTF-8")
  don2<-don[,c(11:44)]
  row.names(don2)<-don$ADM2_PCODE
  
  afc<-CA(don2,ncp = 10,graph = FALSE)
  library(explor)
  res <- explor::prepare_results(afc)
  explor::CA_var_plot(res, xax = 1, yax = 2, lev_sup = FALSE, var_sup = FALSE,
      var_sup_choice = , var_hide = "Row", var_lab_min_contrib = 0, col_var = "Position",
      symbol_var = NULL, size_var = NULL, size_range = c(10, 300), labels_size = 10,
      point_size = 56, transitions = TRUE, labels_positions = "auto", xlim = c(-0.234,
          0.332), ylim = c(-0.187, 0.38))

On enchaîne par une CAH (sur les axes de l’AFC) :

cah <- HCPC(afc,nb.clust = 5,graph = FALSE)
  plot.HCPC(cah,choice = "tree")

Type de structures par âge et sexe des communes du Benin en 2019 (source : OCHA)

plot.HCPC(cah,choice = "map", label.cex=0.3)

Position des classes de la CAH sur les axes de l’AFC (source : OCHA)

Il reste à visualiser les classes. On crée pour cela une fonction assez complexe de comparaison entre le profil de celles-ci et le profil du Bénin, inspirée d’un programme trouvé sur un forum du CIRAD

Profils moyens des classes (source : OCHA)

• la classe 1 : correspond à une population beaucoup plus jeune que la moyenne du pays avec un excédent de 0-9 ans et un déficit de 30 ans et plus. C’est la partie du pays où la transition démographique semble être la plus tardive
• la classe 2 : présente également un profil de population jeune avec un excédent de 5-14 ans, mais se caractérise par une réduction relative des 0-4 ans qui témoigne soit d’une baisse de la fécondité, soit d’un exode rural des jeunes ménages avec enfants.
• la classe 3 correspond au profil moyen du pays.
• la classe 4 est caractérisée par un excédent de personnes âgées et un déficit à la fois de jeunes enfants et de jeunes adultes. Elle correspond à des zones de vieillissment relatif de la population, soit sous l’effet de l’exode rural, soit en raison d’une baisse de la fécondité.
• la classe 5 est caractérise par une forte surreprésentation des jeunes adultes et des persones âgées, associée à un net déficit des enfants jeunes de 0 à 14 ans. Il s’agit vraisemblablement de zones urbaines caractérisées par une fécodité plus tardive et une meilleure espérance de vie.

On peut visualiser la distribution de ces cinq classes sur une carte :

Distribution spatiale des classes (source : OCHA)

2.4.1 Métadonnées

Ce fichier qui est accessible sur HDX en suivant ce lien a été mis à disposition par l’organisation Translators without borders mais l’analyse des métadonnées montre qu’il s’agit en fait de données extraites du recensement géénral de population du Bénin de 2013 :

Created by Translators without Borders, Uploaded on Jun-21, Version 1

Notes and caveats

• All data is drawn from government survey results and is subject to any associated limitations or distortions present in the source data.
• Literacy was measured as the ability to read and write in any language. Languages with population shares less than 0.01% or under 1,000 people (whichever is greater) have been aggregated into the “Other” field.
• All decimal values have been rounded to a maximum of 3 decimal places. As a result language shares may not total 100%.
• Empty values represent non-existent data and should not be treated as zero values. Data is available under an Attribution NonCommercial ShareAlike 4.0 International license (CC BY NC SA 4.0)

Copyright and terms of use

• You are free to share and adapt the data subject to requirements for attribution and non commercial use.
• Any derivative work must be distributed under the same license as the original. Full terms at https://creativecommons.org/licenses/by nc sa/4.0/

Les fichiers fournies comportaient de petites erreurs ou des difficultés de formatage rendant compliquée leur importation dans R. Il a donc fallu les corriger un peu avant de les utiliser. Une fois cela effectué, on a vérifié que les données pouvaient se cartographie en effectuant une jointure avec le fichier administratif au niveau communal.

N.B. nous avons retrouvé l’erreur dans le code de la commune d’Adjara qui est noté BJ0101 alors que le code exact est BJ1001. Nous avons donc reconstruit les fichiers en les harmonisant avec les fonds de carte.

2.4.2 Exemple d’application

Le fichier comporte un ensemble de variables concernant le niveau d’alphabétisation des hommes et des femmes de plus de 6 ans en 2013. Voici à titre indicatif les valeurs par département et pour l’ensemble du pays :

dep<-read.table("data/langues/ben_lang_adm1_correct.csv", sep=";",header =TRUE, dec=",")
  dep<-dep[,c(1:2,77:82)]
  kable(dep)

Ce tableau se prêtera particulièrement bien à des exercices d’initiation à la statistiques et la cartographie sous R puisqu’il comporte à la fois des stocks (population par sexe), des taux (% de population alphabétisée par sexe) et une possibilité de les combiner pour retrouver le nombre de personnes alphabétisée. Il s’inscrit par ailleurs pleinement dans la thématique des inégalités territoriales.

On peut à titre d’exemple construire une cartographie du niveau d’alphabétisation des hommes et des femmes par département du Bénin en 2013 (carte choroplèthe) et superposer par dessus les effectifs correspondants d’hommes et de femmes alphabétisés.

Alphabétisation des hommes et femmes du Bénin par département en 2013 (source : RP Bénin 2013, via OCHA)

On voit apparaître alors une assez forte corrrélation entre les deux distributions, même si les femmes sont moins alphabétisés que les hommes. Ce qui suggère un exercice d’apprentissage de la régression linéaire entre les deux variables.

mod.reg <- lm(dep$literacy_female~dep$literacy_male)
  #stargazer(mod.reg,type = "html")
  ggplot(data = dep) + aes(x=100*literacy_male, y = 100*literacy_female, label = ADM1_NAME) +
                       geom_point(col="red") + 
                       geom_smooth(method="lm") +
                       geom_text(cex=3, nudge_y=2) +
                       scale_x_continuous("Taux d'alphabétisation des hommes (%)") +
                       scale_y_continuous("Taux d'alphabétisation des femmes (%)") +
                       ggtitle("Inégalités d'alphabétisation des départements du Benin par sexe",
                               subtitle = "Source : RP 2013 via OCHA")

Les résultats mettent en évidence le fait que l’alphabétisation des femmes au Bénin en 2013 est en moyenne inférieure de 20% à celle des hommes (le paramètre b de la régression est non significatif) et la relation est statistiquement très significative (r2 = 93 %, p < 0.001). Il existe notamment des résidus c’est-à-dire des départements où la scolarisation des femmes estimée en fonction de celle des hommes est plus forte que prévue (Alibori, Borgou, Littoral) ou plus faible que prévue (Plateau, Couffo, Mono).

Le même modèle appliqué au niveau des communes permettrait une analyse plus fine des résultats à travers une carte des résidus. On pourrait aussi imaginer une analyse mutiscalaire croisant les écarts des communes à la moyenne du pays, de leur département et des communes voisines à l’aide du package mta. Bref, il s’agit sans nul doute d’un excellent exemple pédagogique pour l’école d’été CIST2022.

2.5.1 Métadonnées

Ce fichier qui est accessible sur HDX en suivant ce lien a été mis à disposition par World Food Programme. Cette base qui semble avoir été interrompue présente l’avantage d’offrir une grande profondeur historique et une résolution spatialetrès précise puisque les prix sont localisés par marché et par semaine.

This no longer updated dataset contains Global Food Prices data from the World Food Programme covering foods such as maize, rice, beans, fish, and sugar for 76 countries and some 1,500 markets. It is updated weekly but contains to a large extent monthly data. The data goes back as far as 1992 for a few countries, although many countries started reporting from 2003 or thereafter.

Le fichier contenant les coordonnées de latitude et de longitude des marchés, il peut théoriquement être cartographié facilement à l’aide d’un outil de cartographie dynamique tel que mapview. Mais il peut aussi faire l’objet de multiples agrégations spatiales, temporelles, par produit. Reste évidemment à en apprécier la qualité …

Noter qu’il est nécessaire d’éliminer la seconde ligne du tableau initial sous excel pour pouvoir ensuite charger correctement le fichier. Il est également préférable de lire le fichier avec la fonction fread() du package data.table faut de quoi on peut avoir des problèmes d’encodage plus ou moins délicats à résoudre.

dt<-fread("data/foodprices/wfp_food_prices_ben.csv")
  kable(head(dt))

On peut assez facilement spatialiser le fichier qui comporte en tout 51 marchés assez bien répartiis sur l’ensemble du pays comme le montre la carte ci-desous :

loc<-dt %>% filter(duplicated(market)==FALSE) %>% select(c(2,3,4,5,6)) %>% arrange (admin1,admin2, market)
  markets <- st_as_sf(loc, coords = c("longitude","latitude"))
  st_crs(markets)<- 4326
  #mapview(markets)
  
  
  adm0<-st_read("data/admin/ben_adm0.shp",quiet = TRUE)
  adm1<-st_read("data/admin/ben_adm1.shp",quiet = TRUE)
  adm2<-st_read("data/admin/ben_adm2.shp", quiet = TRUE)
  
  par(mar=c(0,0,0,0))
  plot(adm2$geometry,col="gray90", border = "gray20",lwd=0.5)
  plot(adm1$geometry,col=NA, border = "black",lwd=1, add=TRUE)
  plot(adm0$geometry,col=NA, border = "black",lwd=2, add=TRUE)
  plot(markets$geometry,col="red", bg="yellow",pch=21,add=TRUE,cex=0.7)

2.5.2 Exemple d’application

Cette base de données multidimensionelle peut servir d’exemple pour de nombreux exercices d’apprentissage de la statistique et de la cartographie sous R. Le prix des produits alimnentaires dépend en effet :

• du type de produit (QUOI) ?
• de la date d’achat (QUAND) ?
• du lieu d’achat (OU ?)

Pour chacune des dimensions, on peut réaliser des procédures d’agrégation des valeurs observées ou d’estimation des valeurs manquantes. A titre d’exemple, on extrait de la base de donnée l’ensemble des lignes qui concernent le riz

riz<- dt %>% filter(substr(commodity,1,4)=="Rice") %>% 
                 select(when = date, 
                        where1 = admin1,
                        where2 = admin2,
                        where3 = market,
                        what = commodity,
                        cost = usdprice )
  kable(head(riz))

On peut calculer le prix median du riz sur l’ensemble des points de sondage disponibles en fonction de la catégorie de riz :

tab<- riz %>% group_by(when, what) %>% summarize(median_cost=median(cost), nbobs = n())
  ggplot(tab) + aes(x=when,y=median_cost, color=what, size=nbobs) + geom_point() + 
    scale_x_continuous("Données de prix mensuelles") + 
    scale_y_continuous("Prix médian d'un kilog de riz (en $)")+
    ggtitle(label = "Estimation des prix du riz au Bénin",subtitle = "Source : World Food Program (via HDX)")

Le graphique montre plusieurs choses :

• les données les plus complètes concernent le riz importé qui augmente fortement au moment de la crise de 2007-2008
• la spatialisation des données n’est pas vraiment possible avant 2015 voire 2018
• il y a de fortes différences de prix selon le type de riz

Effectuons maintenant un zoom sur le priz médian du riz importé en 2020-2021 par régions :

tab <- riz %>% filter(when >=as.Date("2020-04-01"), what =="Rice (imported)") %>%
                 group_by(when,where1) %>%
                 summarize(median_cost=median(cost))
  ggplot(tab) + aes(x=when,y=where1, fill=median_cost) +geom_tile() + 
    scale_fill_gradient(low="white", high = "red")+
      scale_x_date("Médiane des valeurs observées par mois") + 
    scale_y_discrete("Département")+
    ggtitle(label = "Prix du riz importé sur les marchés du Benin ($/kg)",subtitle = "Source : World Food Program (via HDX)")

2.6.1 Métadonnées

Ce fichier accessible en cliquant ici concerne le peuplement, c’est-à-dire apparemment l’inventaire de toutes les localisations avec leur nom et leur position en latitude longitude. Il comporte 6306 entrées. Les fichiers sont datés de 2015.

loc<-st_read("data/settlement/ben_plp_NGA.shp", quiet = TRUE)
  par(mar=c(0,0,0,0))
  plot(adm2$geometry,col="lightyellow", border = "gray20",lwd=0.5)
  plot(loc$geometry,col="red",pch=16,cex=0.2, add=T)

Village et peuplement au Benin vers 2015 (Source : OCHA-ROWCA)

2.6.2 Exemple d’application

A titre de vérification de la précision, nous effectuons une superposition sur le fonds de carte OpenStreetMap pour la commune de Ouidah (code BJ0304 ou BEN003004) ou aura lieu l’école d’été du CIST. En cliquant sur les points onpeut comparer leur nom avec celui des localités indiquées par OSM.

map<-adm2[adm2$admin2Pcod=="BJ0304",]
  map<-st_transform(map,4326)
  
  vil<-loc[loc$RowcaCode2=="BEN003004",]
  vil<-st_transform(vil,4326)
  
  
  map <- leaflet() %>% 
              addTiles() %>%
              setView(lat = 6.4, lng=2.1, zoom = 11) %>%
              addPolygons(data = map,
                          fill = FALSE,
                          color = "red",
                          weight = 2) %>%
               addMarkers(data = vil,
                          label = ~FULL_NAME_)
  
  map

Il y a à l’évidence des décalages … Et le contour de la commune lui-même ne semble pas coller exactement avec celui fourni par OSM. Il faudra vérifier laquelle des deux sources est erronée (à moins que ce ne soient les deux ?). Toutefois, la localisation semble rester approximativement juste …

3.1.1 Le tableau “countries”

Il fournit des renseignements de base sur les différents pays, leurs codes, etc.

str(cat$countries)

tibble [304 x 18] (S3: tbl_df/tbl/data.frame)
   $ iso3c             : chr [1:304] "ABW" "AFG" "AFR" "AGO" ...
   $ iso2c             : chr [1:304] "AW" "AF" "A9" "AO" ...
   $ country           : chr [1:304] "Aruba" "Afghanistan" "Africa" "Angola" ...
   $ capital_city      : chr [1:304] "Oranjestad" "Kabul" NA "Luanda" ...
   $ longitude         : num [1:304] -70 69.2 NA 13.2 19.8 ...
   $ latitude          : num [1:304] 12.52 34.52 NA -8.81 41.33 ...
   $ region_iso3c      : chr [1:304] "LCN" "SAS" NA "SSF" ...
   $ region_iso2c      : chr [1:304] "ZJ" "8S" NA "ZG" ...
   $ region            : chr [1:304] "Latin America & Caribbean" "South Asia" "Aggregates" "Sub-Saharan Africa" ...
   $ admin_region_iso3c: chr [1:304] NA "SAS" NA "SSA" ...
   $ admin_region_iso2c: chr [1:304] NA "8S" NA "ZF" ...
   $ admin_region      : chr [1:304] NA "South Asia" NA "Sub-Saharan Africa (excluding high income)" ...
   $ income_level_iso3c: chr [1:304] "HIC" "LIC" NA "LMC" ...
   $ income_level_iso2c: chr [1:304] "XD" "XM" NA "XN" ...
   $ income_level      : chr [1:304] "High income" "Low income" "Aggregates" "Lower middle income" ...
   $ lending_type_iso3c: chr [1:304] "LNX" "IDX" NA "IBD" ...
   $ lending_type_iso2c: chr [1:304] "XX" "XI" NA "XF" ...
   $ lending_type      : chr [1:304] "Not classified" "IDA" "Aggregates" "IBRD" ...

Le tableau comporte 304 observation et il mélange des pays (France), des fragments de pays (Réunion) et des agrégats de pays (Europe). Il faudra donc bien faire attention lors de l’extraction à réfléchir à ce que l’on souhaite utiliser. Par exemple, si l’on veut juste les pays :

## Programme en langage R_base
  # pays<-cat$countries[cat$countries$income_level!="Aggregates",c("iso3c", "country","capital_city","longitude","latitude", "region","income_level")]
  
  
  ## Programme en langage dplyr
  
  pays <- cat$countries %>% 
            filter(income_level !="Aggregates") %>%
            select(iso3c,country, capital_city, latitude, longitude, region, income_level)
  
  
  kable(head(pays))

3.1.2 Le tableau indicators

Il comporte pas loin de 17000 variables … Autant dire qu’il est difficile de l’explorer facilement si l’on ne sait pas ce que l’on cherche.

indic<-cat$indicators
  dim(indic)

[1] 16649     8

kable(head(indic))

vars <- wbsearch(pattern = "CO2",fields="indicator")
  kable(head(vars))

# Programme R-base
  meta<-cat$indicators[cat$indicators$indicator_id %in% c("SP.POP.TOTL","NY.GDP.MKTP.CD","EN.ATM.CO2E.KT"),]
  
  # Programme dplyr
  meta<-cat$indicators %>%
          filter(indicator_id %in% c("SP.POP.TOTL","NY.GDP.MKTP.CD","EN.ATM.CO2E.KT"))
  
  kable(meta)

3.1.3 L’extraction des données

Elle se fait à l’aide de la fonction wb_data qui comporte de nombreuses options.

df   <- wb_data(indicator  = c("NY.GDP.MKTP.CD","SP.POP.TOTL", "EN.ATM.CO2E.KT"))

3.1.4 le choix d’une période de temps

df   <- wb_data(indicator  = c("NY.GDP.MKTP.CD","SP.POP.TOTL", "EN.ATM.CO2E.KT"),
                  start_date = 2017,
                  end_date = 2019)
  dim(df)

[1] 651   7

kable(head(df,6))

df   <- wb_data(indicator  = c("NY.GDP.MKTP.CD","SP.POP.TOTL","EN.ATM.CO2E.KT"),
                  mrv = 1)
  dim(df)

[1] 434   7

kable(head(df,6))

3.1.5 Le choix des unités géographiques

Le paramètre country = permet de choisir les entités spatiales à collecter, soit sous forme de liste de codes, soit à l’aide de valeurs spéciales. Par défaut; il renvoie la liste de tous les pays, mais on peut se limiter à quelques uns seulement à l’aide de leur nom en anglais (risqué …) ou de leur code ISO3 (plus sûr)

df   <- wb_data(indicator  = c("NY.GDP.MKTP.CD","SP.POP.TOTL","EN.ATM.CO2E.KT"),
                  start_date = 2018,
                  end_date = 2018,
                  country = c("BEN","TGO"))
  df$GDP.per.capita <- round(df$NY.GDP.MKTP.CD / df$SP.POP.TOTL,0)
  df$CO2.per.capita <- round(1000*df$EN.ATM.CO2E.KT / df$SP.POP.TOTL,2)
  kable(head(df,6))

df   <- wb_data(indicator  = c("NY.GDP.MKTP.CD","SP.POP.TOTL","EN.ATM.CO2E.KT"),
                  start_date = 2018,
                  end_date = 2018,
                  country = "regions_only")
  df$GDP.per.capita <- round(df$NY.GDP.MKTP.CD / df$SP.POP.TOTL,0)
  df$CO2.per.capita <- round(1000*df$EN.ATM.CO2E.KT / df$SP.POP.TOTL,2)
  kable(df)

3.1.6 Le format de sortie du tableau

Il existe deux façons d’extraire un tableau comprenant plusieurs variables ou plusieurs dates, selon que l’on veut un tableau large (wide) ou étroit. On peut régler la sortie à l’aide du paramètre return_wide qui est TRUE par défaut mais que l’on peut régler sur FALSE.

df   <- wb_data(indicator  = c("NY.GDP.MKTP.CD","SP.POP.TOTL","EN.ATM.CO2E.KT"),
                  return_wide = TRUE,
                  start_date = 2017,
                  end_date = 2018,
                  country = c("BEN","TGO"))
  df

# A tibble: 4 x 7
    iso2c iso3c country  date EN.ATM.CO2E.KT NY.GDP.MKTP.CD SP.POP.TOTL
    <chr> <chr> <chr>   <dbl>          <dbl>          <dbl>       <dbl>
  1 BJ    BEN   Benin    2017           7360   12701654743.    11175192
  2 BJ    BEN   Benin    2018           7910   14262407012.    11485035
  3 TG    TGO   Togo     2017           2070    6395472574.     7698476
  4 TG    TGO   Togo     2018           2260    7112200725.     7889095

df   <- wb_data(indicator  = c("NY.GDP.MKTP.CD","SP.POP.TOTL","EN.ATM.CO2E.KT"),
                  return_wide = FALSE,
                  start_date = 2017,
                  end_date = 2018,
                  country = c("BEN","TGO"))
  df[,1:7]

# A tibble: 12 x 7
     indicator_id   indicator          iso2c iso3c country  date        value
     <chr>          <chr>              <chr> <chr> <chr>   <dbl>        <dbl>
   1 NY.GDP.MKTP.CD GDP (current US$)  BJ    BEN   Benin    2018 14262407012.
   2 NY.GDP.MKTP.CD GDP (current US$)  BJ    BEN   Benin    2017 12701654743.
   3 NY.GDP.MKTP.CD GDP (current US$)  TG    TGO   Togo     2018  7112200725.
   4 NY.GDP.MKTP.CD GDP (current US$)  TG    TGO   Togo     2017  6395472574.
   5 SP.POP.TOTL    Population, total  BJ    BEN   Benin    2018    11485035 
   6 SP.POP.TOTL    Population, total  BJ    BEN   Benin    2017    11175192 
   7 SP.POP.TOTL    Population, total  TG    TGO   Togo     2018     7889095 
   8 SP.POP.TOTL    Population, total  TG    TGO   Togo     2017     7698476 
   9 EN.ATM.CO2E.KT CO2 emissions (kt) BJ    BEN   Benin    2018        7910 
  10 EN.ATM.CO2E.KT CO2 emissions (kt) BJ    BEN   Benin    2017        7360 
  11 EN.ATM.CO2E.KT CO2 emissions (kt) TG    TGO   Togo     2018        2260 
  12 EN.ATM.CO2E.KT CO2 emissions (kt) TG    TGO   Togo     2017        2070 

3.2.1 Natural Earth

Nous allons ici utiliser le fonds de carte Natural Earth qui est un fonds de carte libre de droit et mis à jour régulièrement. Le site web du projet se situe à l’adresse suivante :

https://www.naturalearthdata.com/

Il indique ses objectifs comme suit :

“Natural Earth is a public domain map dataset available at 1:10m, 1:50m, and 1:110 million scales. Featuring tightly integrated vector and raster data, with Natural Earth you can make a variety of visually pleasing, well-crafted maps with cartography or GIS software.[…] Natural Earth was built through a collaboration of many volunteers and is supported by NACIS (North American Cartographic Information Society), and is free for use in any type of project (see our Terms of Use page for more information).”

On peut télécharger les différents fonds de carte sur le site web, mais dans une perspective de mise à jour automatique régulière du fonds de carte il est plus pertinent d’utiliser l’API rnaturalearthqui permet d’accéder directement à la plupart des fonds de carte avec juste quelques lignes de code. Il suffit pour cela de commencer par installer et charger le package.

library("rnaturalearth")
  library("rnaturalearthdata")

3.2.2 le fonds de carte countries110 (175 unités)

On va télécharger tout d’abord le fonds de carte des pays du Monde avec une forte généralisation des contours countries110 et le transformer en objet de type spatial feature du package sf du package avant de le visualiser et d’ examiner le nombre d’unités

map<-st_as_sf(countries110)
  class(map)

[1] "sf"         "data.frame"

ggplot(data = map) +
          geom_sf(fill="lightyellow") +
          theme_bw()

Ce fonds de carte comporte 175 unités spatiales, mais de quoi s’agit-il exactement. Les métadonnées associées permettent de se faire une idée plus précise de la nature exacte de ces unités. Prenons pour cela quelques exempes

sel<-map[map$adm0_a3 %in% c("FRA", "NCL","ATA","ATF","USA", "PRI","CHN","TWN","MAR", "SAH","CHN","TWN","ISR","PSX"),c("sovereignt","sov_a3","type","admin", "adm0_a3","name","note_adm0","iso_a3","wb_a3")]
  kable(sel)

Les exemples présentés dans le tableau ci-dessus montrent la complexité du problème de définition et de représentation cartographique des “pays” ou “bouts du monde”. Quelques remarques :

1. La France (FR1) en tant qu’état souverain regroupe ici cartographiquement la partie métropolitaine du pays et les Départements d’Outre-Mer (Guyane Française, Réunion, Martinique, Guadeloupe) en une seule entité spatiale, mais elle met à part la Nouvelle Calédonie et les îles antarctiques.
2. Porto Rico (PRI) est considéré comme une dépendance des Etats-Unis (US1) au même titre que la Nouvelle Calédonie(NCL) est considérée comme une dépendance de la France (FR1).
3. Le Sahara occidental (SAH) est considéré comme une zone indéterminée bien qu’il soit occupé par le Maroc (MAR).
4. la Palestine (PSX) est considéré comme une zone disputée mais rattachée en terme de souveraineté à Israël (ISR) et une note précise qu’elle est partiellement semi-administrée. Le code sur trois caractères des territoires palestiniens est très variable selon les organisations (PSX, PSE, WBG).
5. Taïwan (TWN) est présenté comme un état souverain, mais son code ISO3 est manquant pour la banque mondiale car la Chine refuse de le reconnaître.
6. Plusieurs états souverains de petite taille sont absents de ce fonds de carte qui ne regroupe que 175 unités soit moins que les 193 pays membres des Nations-Unies. La plupart des îles du Pacifique sont en particulier éliminées car leur surface les rendrait invisible pour le degré de généralisation cartographique adopté.

3.2.3 le fonds de carte sovereignty110 (171 unités)

On peut obtenir un fonds différent en installant le package complémentaire rnaturalearthdata qui permet notamment de distinguer le fonds de carte des countries (c’est-à-dire des “bouts du monde” souverains ou non) et des sovereignty (c’est-à-dire des états souverains)

library(rnaturalearthdata)
  map<-st_as_sf(sovereignty110)
  
  ggplot(data = map) +
   geom_sf(fill="lightyellow") +  
    theme_bw()

Le fonds de carte permet désormais de récupérer la plupart des pays souverains du Monde, y compris les petits états insulaires du Pacifique, mais il fait disparaître de façon sélective les territoires indéterminés ou disputés. Ainsi, le Sahra Occidental demeure partiellement séparé du Maroc mais les territoires palestiniens sont annexés à Israël ainsi que le plateau du Golan ce qui n’est évidemment pas un choix neutred’un point de vue géoolitique.

par(mfrow=c(1,2))
  
  map2<-map %>% filter(sov_a3 %in% c("ISR","JOR","SYR","LBN","EGY"))
  ggplot(data = map2) +
      geom_sf(fill=c("gray80","orange","gray80","gray80","gray80")) +
      ggtitle("Limits of Israël") +
    theme_minimal()

map3<-map %>% filter(sov_a3 %in%c("MAR","SAH","DZA","MRT"))
  ggplot(data = map3) +
      geom_sf(fill=c("gray70","orange","gray70","lightyellow")) +
      ggtitle("Limits of Morocco") +
    theme_minimal()

3.2.4 Le fonds de carte countries50

On peut également choisir un fonds moins généralisé dans lequel tous les petits pays seront présents

map<-st_as_sf(countries50)
  ggplot(data = map) +
      geom_sf(fill="lightyellow") +
    theme_bw()

Il existe toute une série d’autres fonds de carte dans le package Natural Earth, notamment avec des résolutions plus précises, mais on se limitera ici à l’exploration des fonds de carte utile pour produire des cartes à contour généralisé couvrant le monde entier.

3.3.1 Etape 1 : récupération des données statistiques

Nous commençons par récupérer les données brutes de population et de CO2 en 2018 pour les pays de la CEDEAO et on y ajoute les latitues et longitues des centres des pays.

cedeao<-c("BEN","BFA","CPV","CIV", "GMB","GIN","GNB","GHA","LBR","MLI","NER","NGA","SEN","SLE","TGO")
  df   <- wb_data(indicator  = c("SP.POP.TOTL", "EN.ATM.CO2E.KT"),
                  return_wide = TRUE,
                  start_date = 2018,
                  end_date = 2018,
                  country = cedeao)
  kable(df)

Nous renommons les variables pour avoir un tableau plus simple ou la population est en millions d’habitants, les émissions de CO2 en millions de tonnes. On y ajoute l’intensité des émissions en tonnes par habitant.

don <-df %>% select(ISO3 = iso3c, NOM = country, POP = SP.POP.TOTL, CO2 = EN.ATM.CO2E.KT) %>%
              mutate(POP = POP/1000000, CO2 = CO2/1000, CO2_hab = CO2/POP)
  kable(don,digits = 2)

3.3.2 Etape 2 : Récupération du fonds de carte

on récupère ensuite le fonds de carte en ne gardant que les pays de la CEDEAO

cedeao<-c("BEN","BFA","CPV","CIV", "GMB","GIN","GNB","GHA","LBR","MLI","NER","NGA","SEN","SLE","TGO")
  map<-st_as_sf(countries110)
  map<-st_as_sf(countries110) %>% 
          select(adm0_a3,geometry) %>% 
          rename(ISO3 = adm0_a3) %>%
          filter(ISO3 %in% cedeao)
  plot(map$geometry, col="lightyellow")

3.3.3 Etape 3 : Jointure du fonds de carte et des statistiques

mapdon <- right_join(don,map) %>% st_as_sf()
  kable(head(mapdon))

3.3.4 Etape 4 : Visualisation avec mapsf

mf_theme("agolalight")
  mapdon %>% 
    mf_map("CO2_hab", 
           "choro",
           breaks="jenks",
           leg_pos="bottomleft",
           leg_title = "en tonnes/hab.")%>%
    mf_map("CO2", 
           "prop",
           col="red",
           leg_pos = "topleft",
           leg_title = "en millions de tonnes")
  mf_title("Emissions de CO2 des pays de la CEDEAO en 2018")